/*
 *  VectorUtils.h
 *  NSE_SOLVER
 *
 *  Created by Patrick Zulian on 1/14/11.
 *  Copyright 2011 ICS - USI - Lugano. All rights reserved.
 *
 */

#ifndef NSE_SOLVER_VECTOR_UTILS_H
#define NSE_SOLVER_VECTOR_UTILS_H

#include "Base.h"


namespace nsesolver {
	
	
	template< class VectorType >
	const VectorType Absolute(const VectorType &v)
	{
		using std::abs;
		
		VectorType result;
		for (int i = 0; i < v.size(); i++) {
			result[i] = abs(v[i]);
		}
		
		return result;
	}
	
	template< class VectorType >
	const VectorType CanonicalBasisVector(const int spanningDim)
	{
		VectorType result;
		result.allSet(0);
		result[spanningDim] = 1;
		return result;
	}
	
	
	template< class Vector3Type >
	const Vector3Type Cross(const Vector3Type &v1, const Vector3Type &v2)
	{
		NSE_SOLVER_ASSERT( v1.size() == 3 );
		NSE_SOLVER_ASSERT( v2.size() == 3 );
		
		
		Vector3Type result;
		result[X] = (v1[Y] * v2[Z]) - ( v2[Y] * v1[Z] );
		result[Y] = (v1[Z] * v2[X]) - ( v2[Z] * v1[X] ); 
		result[Z] = (v1[X] * v2[Y]) - ( v2[X] * v1[Y] );
		
		return result;
	}
	
	//template< class Vector4Type >
//	const Vector4Type Cross(const Vector4Type &v1, const Vector4Type &v2, const Vector4Type &v3)
//	{
//		Vector4Type result;
//		result[0] = 
//	}
	
	
	
	template< class VectorType, typename T >
	const VectorType MeshVec1(const T &h0, const T &by, const T &hEnd)
	{
		int n = ceil( (hEnd - h0)/(by) + 1.0 );
		VectorType result(n);
		
		T h = h0;
		for (int index = 0; index < n; index++) {
			result[index]= h;
			h += by;
		}
		return result;
	}
	
	template< class VectorType, typename T >
	const VectorType MeshVec1(T h0, T hEnd)
	{
		return MeshVec1< VectorType, T >(h0, 1, hEnd);
	}
		
}

#endif
